%--------------------------------------------------------------------------
% computes the eigenvalues of the adjoint problem using a 'frozen time'
% approach
%
% assumes the basic state is in the t = O(1 / delta) regime
%--------------------------------------------------------------------------


function M = comp_adj(k, p, h0)

if nargin == 1
    p = params;
end

if nargin < 3
    h0 = 1;
end

% basic state at z = h_bar
v_0 = (1 - p.beta) * (1 - 1 / h0);
c_b = p.beta + v_0 - 1/2 * p.delta * (1 - p.beta) * (p.beta + v_0);

N = p.N;
z = linspace(0, h0, N);
h = z(2) - z(1);

M = spalloc(N, N, 4 * N);


        
if (k < 0.01)
    w = @(z) p.Ma * z.^2 .* (k^2 / 4 * (1 - z) + p.delta / p.Ma / 2 * (3 - z));
else
    A = w_coeffs(1, h0, k, p);
    w = @(z) (A(1) * z + A(2)) .* cosh(k*z) + (A(3) * z + A(4)) .* sinh(k*z);
end
c1bz = @(z) -p.delta * (1 - p.beta) * (p.beta + v_0) * z / h0^2;

% build the matrix

for i = 1:N
    
    if(i == 1)
        M(i,i) = -2 / h^2 - k^2;
        M(i,i+1) = 2 / h^2;
    elseif (i == N)
        M(i,1) = -w(z(1)) * c1bz(z(1));
        M(i,2:N-2) = -2 * w(z(2:N-2)) .* c1bz(z(2:N-2));
        M(i,N-1) = 2 / h^2 - 2 * w(z(N-1)) * c1bz(z(N-1));
        M(i,N) = -2 / h^2 - 2 * p.delta * (1 - 2*c_b) / h - w(z(N)) * c1bz(z(N)) - k^2;
    else
        M(i,i) = -2 / h^2 - k^2;
        M(i,i-1) = 1 / h^2;
        M(i,i+1) = 1 / h^2;
    end
end

